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Abstract. 

We use a time-dependent dynamical mean-field-hydrodynamic model to study the 
formation of fermionic dark solitons in a trapped degenerate fermi gas mixed with a 
Bose-Einstein condensate in a harmonic as well as a periodic optical-lattice potential. 
The dark soliton with a "notch" in the probability density with a zero at the minimum 
is simulated numerically as a nonlinear continuation of the first vibrational excitation 
of the linear mean-field-hydrodynamic equations, as suggested recently for pure bosons. 
We study the free expansion of these dark solitons as well as the consequent increase in 
the size of their central notch and discuss the possibility of experimental observation 
of the notch after free expansion. 

PACS numbers: 03.75.Lm, 03.75.Ss 
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1. Introduction 

Due to a strong repulsive Pauli-blocking interaction at low energies among spin-polarized 
fermions, there cannot be an evaporative cooling leading to a quantum degenerate 
fermi gas (DFG) [1]. Trapped DFG has been achieved only by sympathetic cooling 
in the presence of a second boson or fermion component. Recently, there have been 
successful observation [1-4] and associated experimental [5-7] and theoretical [8-13] 
studies of degenerate boson-fermion mixtures by different experimental groups [1-4] 
in the following systems: 6 ' 7 Li [3], 23 Na- 6 Li [4] and 87 Rb- 40 K [5,6]. Also, there have 
been studies of a degenerate mixture of two components of fermionic 40 K [1] and 6 Li [2] 
atoms. The collapse of the DFG in a boson-fermion mixture 87 Rb- 40 K has been observed 
and studied by Modugno et al. [5,10,13]. In these studies of a mixed Bose-Einstein 
condensate (BEC) and a DFG the initial states were the stationary ground states of the 
systems. 

In this paper we study the possibility of the formation of fermionic dark solitons in a 
mixture of a DFG with a BEC using a coupled time-dependent mean-field-hydrodynamic 
model where the bosonic component is treated by the mean-field Gross-Pitaevskii (GP) 
equation [14] and the fermionic component is treated by a hydrodynamic model [11,12]. 
This time-dependent mean-field-hydrodynamic model was suggested recently by the 
present author [13] to study the collapse dynamics of a DFG and is a time-dependent 
extension of a time-independent model suggested for the stationary states by Capuzzi 
et al. [11,12] based essentially on a Thomas- Fermi- Weizsacker approximation. 

Zakharov and Shabat [15] have shown that the dimensionless nonlinear Schrodinger 
(NLS) equation in the repulsive or self-defocusing case [16] 

iu t + u xx - \u\ 2 u = 0. (1.1) 

sustains the following dark and grey solitons [17]: 

u(x, t) = r(x — ct) exp[— i{<p(x — ci) — fit}], (1.2) 

with 

r 2 (x — ct) = 1] — 2/t 2 sech 2 [/t(x — ct)], (1.3) 

(f>(x — ct) = tan _1 [— 2k/c tanh{/t(x — ct)}], (1-4) 

k = vW- c 2 )/2, (1.5) 

where c is the velocity of the soliton, fi the parametric energy, and rj related to intensity. 
Soliton (jl.2j) having a "notch" over a background density is grey in general. It is dark 
if density \u\ 2 = at the minimum. The soliton can move freely with velocity c and at 
zero velocity the soliton becomes a dark soliton: \u(x, t)\ = y^rj tanh[x \Jrj/2]. 

The similarity of the NLS equation (jl.lj) to the GP equation (|2.1|) (below) imply the 
possibility of a dark soliton in a trapped BEC [18]. It has been suggested that the dark 
soliton of a trapped BEC could be a stationary eigenstate of the GP equation [18, 19] 
as in the case of the trap-less NLS equation. The usual search for a dark soliton in 
the GP equation proceeds through time evolution starting with the ansatz [20, 21] 
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M DS = tanh(x)iiQ(x), where uq is the ground state of the GP equation and the tanh(x) 
factor is introduced in analogy with the dark soliton in the NLS equation. However, 
the function wj^g is not an eigenfunction of the GP equation and hence this procedure 
leads to numerical instability on time evolution [20,21]. It has been demonstrated [22] 
that the time evolution of the GP equation with the initial state wpjg leads to a dark 
soliton which is the lowest vibrational excitation of the system. Exploring this, a stable 
numerical procedure has been suggested [22] for the simulation of dark soliton which we 
use in this investigation. 

We study the formation of fermionic dark solitons in a DFG mixed with a BEC in 
a harmonic as well as a periodic optical-lattice trap. There have been experimental [23] 
and theoretical [24] studies of the formation of dark solitons in a harmonically trapped 
BEC. In view of this, here we study for the first time the possibility of the formation of 
a fermionic dark soliton in a DFG using a mean-field model for a mixture of DFG and 
BEC. 

Collective excitations in the form of solitons and vortices in trapped fermions have 
also been investigated recently by Damski et al. [25] and Karpiuk et al. [26]. However, 
they considered isolated ultra-cold fermions and not a realistic mixture of trapped DFG 
and BEC as in the experiments and as discussed in this paper. Also they did not 
demonstrate the existence of stable dark solitons with a zero at the center of the notch 
as in the present study. They identified grey soliton-like structure with a shallow dip in 
the isolated fermionic density distribution quite distinct from the stable fermionic dark 
solitons in a DFG-BEC mixture as noted in this investigation. 

Though the present dark solitons are numerically stable in the mean-field 
formulation, they could be unstable physically due to quantum fluctuations [27]. The 
effect of quantum fluctuations is lost in the mean-free model and can only be studied in a 
field-theoretic approach. Moreover being an excited state they are thermodynamically 
unstable. There have been suggestions about how to excite a dark soliton by phase 
imprinting method [25,26]. The dark soliton is the lowest vibrational excitation of the 
BEC [22] and there have also been investigations about how to attain such excited 
states [28]. Nevertheless, despite different suggestions about how to excite a dark 
soliton [25,26,28], experiments to date have not yet generated a stationary dark soliton. 
Experimentally, so far the dark solitons have been unstable [18, 19]. Although we 
cannot make definite suggestion(s) for the formation of stable fermionic dark solitons, 
considering that they are stationary excitations of the mean-field equation their creation 
seems possible at least as a non-stationary dark soliton which may turn grey and oscillate 
before decaying due to quantum fluctuations and thermodynamic effects. 

Lately, the periodic optical-lattice potential has played an essential role in many 
theoretical and experimental studies of Bose-Einstein condensation, e. g., in the study of 
Josephson oscillation [29] and its disruption [30], interference of matter- wave [31], BEC 
dynamics on periodic trap [32], etc. The periodic optical-lattice confinement generated 
experimentally by a standing-wave laser field creates a BEC in an entirely different shape 
and trapping condition form a conventional harmonic oscillator trapping. In view of this 
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we study the possibility of the formation of a fermionic dark soliton in an optical-lattice 
potential. The formation of a bosonic dark soliton in an optical-lattice potential has 
already been investigated [21]. 

The central notch is the earmark of a dark soliton. Experimentally, a dark soliton is 
identified after removing the traps so that a free expansion of the DFG allows the notch 
to widen and be photographed clearly. In view of this we study a free expansion of dark 
solitons in a DFG-BEC mixture and study the possibility of detection of a fermionic 
dark soliton in the laboratory. 

In section 2 we present an account of the time-dependent mean-field model 
consisting of a set of coupled partial differential equations involving a BEC and a DFG. 
In the case of a cigar-shaped system with stronger radial trapping, the above model is 
reduced to an effective one-dimensional form appropriate for the study of dark solitons. 
In section 3 we present our results for stationary fermionic dark solitons as well as a 
study of their free expansion in a boson-fermion mixture. Finally, a summary of our 
findings is given in section 4. 



2. Nonlinear mean- field- hydrodynamic model 



The time-dependent BEC wave function \l/(r, t) at position r and time t may be described 
by the following mean-field nonlinear GP equation [14] 



ih 



d h 2 v 2 r 



+ Vs(r) + g B Bn B Vb(t, t) = 



(2.1) 



dt 2rriB 

with normalization J dr|^s(r, t)\ 2 = Nb- Here tub is the mass and Nb the number of 
bosonic atoms in the condensate, n# = |^/s(r,t)| 2 is the boson probability density, 
Qbb = 47r^ 2 asB/m j B the strength of inter- atomic interaction, with a B B the boson- 
boson scattering length. The trap potential with axial symmetry may be written as 
V B (r) = \m B oo 2 (p 2 + v 2 z 2 ) where uo and voj are the angular frequencies in the radial (p) 
and axial (z) directions with v the anisotropy parameter. The probability density np 
of an isolated DFG in the Thomas-Fermi approximation is given by [10] 
[max(0,{6 F -W(r)})] 3 /2 



Up 



A 3/2 



(2.2) 



where A = ft 2 (6vr 2 ) 2 / 3 / (2m F ), ep is the Fermi energy, rrtp is the fermionic mass, and the 
function max denotes the larger of the arguments. The confining trap potential Vp(r) 
has axial symmetry. The number of fermionic atoms iVV is given by the normalization 
J drriF = Np. 

We developed a set of practical time-dependent mean-field-hydrodynamic equations 
for the interacting boson-fermion mixture starting from the following Lagrangian density 
[13] 



-h 



~dt 



dt 



dt 



+ 



h 2 \v r m B \ 

2m B 



+ V b \Vb\ 2 + t;9bb\Vb\ 



5/3 
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h 2 \V r Jn~^\ 2 . . 3.. 

— - h Vf njr + -A\n F 

brriF 5 

+ g BF n F \^ B \ 2 , (2.3) 

where qbf = 2irh 2 a B F /thr with the boson-fermion reduced mass m B = m B rriF / '{m B + 
mp)i where a B p is the boson-fermion scattering length. 

It may not be entirely proper to define an average fermionic wave function 

— V™f * n a DFG like in a BEC. The correct fermionic wave function is to be 
calculated from a Slater determinant Schrodinger equation for degenerate fermions [9]. 
However, the probability density rtp of a DFG calculated in this fashion should lead to 
reasonable results [9] and has led to proper probability distribution for a DFG [12, 13] 
as well as results for collapse of a DFG [10, 12, 13] in agreement with experiment. This 
approach has also been used successfully to predict a fermionic bright soliton in a boson- 
fermion mixture [33]. 

The terms in the first round bracket on the right-hand side of ()2.3|) are the standard 
Gross-Pitaevskii terms for the bosons and are related to a Schrodinger-like equation [14]. 
However, terms in the second round bracket, although bear a resemblance with the 
first, are derived from the hydrodynamic equation of motion of the fermions including a 
Weiszacker kinetic energy and are not related to a Schrodinger-like equation [11]. Hence, 
the second kinetic energy term has a different mass factor 6mp and not the conventional 
Schrodinger mass factor 2m B as in the first term. Finally, the last term in this 
equation corresponds to an interaction between bosons and fermions. The interaction 
between bosons and between bosons and fermions are described by contact potentials 
parametrized by coupling constants g BB and g BF defined above. The interaction between 
fermions in spin polarized state is highly suppressed due to Pauli blocking and has been 
neglected in ()2.3|) and will be neglected throughout this paper. 

Recently, Jezek et al. [34] used the Thomas-Fermi- Weizsacker kinetic energy 
term Tp of fermions in their formulation which, in our notation, will correspond 
to a fermionic kinetic energy of % 2 \V ry /np\ 2 / (Qmp) in ()2.3j) in place of the present 
term h 2 \ V r y/np\ 2 / {Qm F ). This kinetic energy term contributes little to this problem 
compared to the dominating 3A\n F \ 5 ^ 3 /5 term in ()2.3j) and is usually neglected in the 
Thomas-Fermi approximation. However, its inclusion leads to a probability density 
which is a smooth and analytic function of the space variable [34]. For a discussion of 
these two fermionic kinetic energy terms we refer the reader to [11,34,35]. 

With the Lagrangian density ()2.3|) . the Euler-Lagrange equations of motion become 

[13]: 



d h 2 v 2 

ih— - - — L + V B (r) + g B Bn B + ^fI^fI 
ot 2m B 



tt fl (r,t) = 0, (2.4) 



d h 2 v 2 1 

-ih—-—^ + VF{Y) + A\n F \ 2 l z + g B Fn B y/n F (r, t) = 0. (2.5) 
ot brriF -I 

When the nonlinearity in ()2.5|) is very large, the kinetic energy term in this equation 
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can be neglected and the time-independent stationary form of this equation becomes 

[max(0, {e F - V F (r) - g B Fn B })] 3/2 fn R . 

n F -p^ , (2.b) 

which is the generalization of (|2.2jl in the presence of boson-fermion coupling. Equation 
()2.6|) has been used by Modugno et al. [10] for an analysis of a BEC coupled to a DFG. 
We shall see in the following that in actual experimental condition the nonlinearity in 
()2.5|) is quite large and ()2.6|) is a good approximation. We note that the Lagrangian 
density of the formulation of Jezek et al. [34] reduces to the present Lagrangian density 
in this approximation upon the neglect of the fermionic kinetic energy term. 

The solution of the coupled three-dimensional equations above for studying dark 
solitons in a boson-fermion mixture is a formidable task. Hence, we shall reduce f!2.4j) 
and ()2.5j) to the minimal one-dimensional form suitable for the study of dark solitons in 
a cigar-shaped geometry for v << 1. We perform this reduction below where we take 
Vg(r) = V F (r) = \rriB<jJ 2 (p 2 + v 2 z 2 ) + Vo sin 2 (27T2;/A) which corresponds to a suppression 
of uj and vuj for fermions by a factor ^JtobJ^f as in the study by Modugno et al. [10] and 
Jezek et al. [34]. In the confining potential we also include the following optical-lattice 
potential: V Q sm 2 {2tx z / X) [36]. Here V is the strength of the optical-lattice potential 
and A is the wave length of the laser. 

For v << 1, (|2.4j) and ()2.5|) can be reduced to an effective one-dimensional 
form by considering solutions of the type ^(r, t) = 4>b{z, t)tp^ (p) and y/n F (r, t) = 
(f) F (z,t)il)P(p), where 

i ,(o)/ N.2 _ / M t up 2 \ 



nh \ h 

corresponds to the respective circularly symmetric ground state wave function in the 
absence of nonlinear interactions and satisfies 

- ^> { b ] + ^ 2 M 0) = (2-8) 

ft 2 2 , (0) . 1 2 2, (0) / mB k , (0) 

-6^ V ^ + 2 mBUJ P ^ = ' 

(2.9) 

with normalization 2tt J °° \%l)f\p)\ 2 pdp = 1. Now the dynamics is carried by <fii(z,t) 
and the radial dependence is frozen in the ground state ijjf\p). The separation of the 
variables is suggested by the structure of (|2.4|) and (|2.5jl . 

Averaging over the radial mode ipf\p), i.e., multiplying (|2.4j) and ()2.5|) by ip^*(p) 
and integrating over p, we obtain the following one-dimensional dynamical equations 
[37]: 

d h 2 d 2 1 9 2 2 2 ( 

— in— — r + -rriBV uj z + Vq sin 



dt 2m B dz 2 2 u V x 

+ F BB \<j) B \ 2 + FBF\M 2 ]Mz,t) = 0, (2.10) 
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h 2 



lh dt 



d 2 1 

i 2 2 2 , -tr ■ 2 

7- + tP-bv oj z + V sin 

bm F oz A 2 



2nz 



+ F FF \<p F \ i/3 + Fbf\M 2 Mz, t) = 



where 



BB 



BF 



' FF 



9bb 



9bf 



2^K' 



SZ°\1>P\ 2 



\^\ 2 pdp 



r 

, Jo 



|</4 0) | 2 + 4 /Vp 



Mbf^ 



h 



|V4° \ 2 pdp 



3A 
~5~ 



M f uj 



2/3 



(2.11) 



(2.12) 



(2.13) 



(2.14) 



Jo l^-B 

In (12.10)) and (J2.ll)) we have included the optical-lattice potential and the normalization 
there is given by £^ \(f>i(z, t)\ 2 dz = N { . 

For calculational purpose it is convenient to reduce the sets (12.10)1 and 1)2.11)1 to 
dimensionless form by introducing convenient dimensionless variables. Although the 
algebra is quite straightforward, the expressions become messy with different factors of 
masses. As we shall not be interested in a particular boson-fermion system in this paper, 
but be concerned with the formation of fermionic dark solitons in general, we take in the 
rest of this paper mg = 3m F = m( 87 Rb), whence ttir = 3m F /4:, Mb = Mp — tub, and 
Mbf — ^s/2 and where m( 87 Rb) is the mass of the Rb atom. In the two experimental 
situations of [4,5] m B ~ 3m F . 

In (J2.10)) and (J2.ll)) . v « 1 and we consider the dimensionless variables r = tvu/2, 
y = z/l z ,Xi = \/{l z /Ni)(f)i, with l z = ^h/(uajm B ), so that 
■ d d 2 2 . 2 (2-ny 

|2 



+ -A BB \Xb\ 



+ N B f \Xf\ 



Xb(v,t) 



0. 



(2.15) 



d 2 
dy 2 

+ N FF | X f| 4/3 



'2vry 
An 



+ -A pb \xb\ 



Xf(v,t) = 0, 



(2.16) 



where Nbb 



(i/p)a BB N B /l z , N BF = (8/v)a BF N F /l z , N FB = (8/u)a BF N B /l z , and 
N FF = 9(6iiN F /h') 2 / 3 /5. In (12.15)) and (j2.16j) . the normalization condition is given by 
JZ° \Xi{Vi T )\ 2 dy = 1 and t>n = 2Vo/(huju) is the reduced strength of the optical-lattice 
potential and A = A// 2 is the dimensionless wave length. 



3. Numerical Result 



The coupled mean-field-hydrodynamic equations ()2.15|) and (J2.16)) for dark solitons 
are solved numerically using a time-iteration method based on the Crank-Nicholson 
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discretization scheme elaborated in [38]. We discretize the mean-field- hydrodynamic 
equations using time step 0.0005 and space step 0.025. 

We performed the time evolution of the set of equations (j2.15j) and ()2.1ti|) setting 
Nbb = N BF = N FB = N FF = vq = and starting with the eigenfunction of the 
lowest excited state of the linear harmonic oscillator problem as suggested recently [22], 
e.g., with Xb(v, t ) = Xf(v,t) = A/2vr~ 1//4 |/exp(— y 2 /2) exp(— 3ir). During the course 
of time evolution the nonlinear terms are switched on very slowly and the resultant 
solutions iterated (about 50000 times) until convergence was obtained. If converged 
solutions are obtained, they correspond to the required dark solitons. In the present 
approach the time evolution starts with and proceeds through successive eigenfunctions 
of the coupled mean- field equations (j2.15|) and f)2. 16j) . Hence it leads to stable numerical 
results. The usual numerical procedure for the calculation of the dark solitons starts 
with an approximate solution of the mean-field equation and hence leads to numerical 
instability on time evolution [20-22]. 

We solve (j2.15j) and (J2.16)) for dark solitons. In this case the nonlinearity N FF 
could be very large for N F > 100, which may require special care for obtaining accurate 
numerical solutions. In our calculation we use v = 0.1, Vq = 0, vu — 27T x 100 Hz and 
rriB to be the mass of 87 Rb. Consequently, / ~ 1 fim and unit of time r = IjivuS) is 3 
ms. We also take N F = 100, N B = 100, l z = 1 //m, and a B B = a BF = 5 nm. With these 
parameters the nonlinearities in ()2.15j) and ()2.16|) are Nbb = 20, Nbf = 40, N F b = 40, 
and N FF = 1275. 

The results of the present study on dark solitons are presented in figures 1. In figure 
1 (a) we plot the stationary functions \<f>(z, t)\ of the boson and fermion dark solitons. 
Both functions have a notch in the middle. The notch in the bosonic function is wider 
than that in the fermionic function. In this case the fermionic nonlinearity (1275) is 
large. Consequently, the extension of the fermionic function in figure 1 (a) is also much 
larger than the bosonic function. 

One way to observe the dark solitons experimentally is to allow them to expand 
freely while the central notch will become wider in size to be visible and photographed. 
With this in mind we study the free expansion of the boson-fermion mixture. The 
snapshots of the bosonic and fermionic functions at regular intervals of time during this 
expansion are shown in figures 1 (b) and (c), respectively. After expansion, the central 
notch in the bosonic dark soliton expands as we find from figure 1 (b) but the notch in 
the fermionic dark soliton in figure 1 (c) does not expand enough to be visible. This 
will make the fermionic dark soliton more difficult to observe experimentally. However, 
this behavior is quite expected in a freely expanding fermionic dark soliton with a large 
nonlinear repulsion. In the fermionic equation both the fermion-fermion and fermion- 
boson interactions are highly repulsive. The dark soliton has a notch (hollow region) 
at the center. Because of the very strong repulsion, the DFG tends to expand in all 
directions including the radially inward direction to fill the central hollow space as 
well as outward directions. This inward repulsive force balances partially the outwards 
kinetic pressure and does not allow the central notch to expand substantially during 
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Figure 1. (a) The stationary function \<fi(z, t) \ for bosonic and fermionic dark solitons 
vs. z with Np = Nb = 100, v = 0.1 and cibb — o-bf = 5 ran. We show in (b) and (c) 
the profiles of the bosonic and fermionic functions |0s(z,t)| and \4>F(z,t)\, respectively, 
of the degenerate mixture (a) during free expansion at regular intervals of time. The 
nonlinearities are N B b = 20, N B f = 40, N FB = 40, Npp — 1275. 
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Figure 2. The profiles of the oscillating function |0(z,t)| for (a) bosonic and (b) 
fermionic dark solitons of figure 1 (a) vs. z. The oscillation was originated by 
jumping the boson-fermion scattering length so that the nonlinearities Nbf and Nfb 
are suddenly jumped from 40 to 44. 



free expansion so as to be easily observable. This is not the case for a moderately 
repulsive or attractive pure single-component BEC, where the outward kinetic pressure 
overcomes any repulsion among the atoms and the BEC expands only in the radially 
outward direction with a widening of the notch. In the bosonic wave function of figure 1 
(b) the notch expands reasonably during the expansion of the BEC. However, during the 
expansion in figures 1 (b) and (c) the central notch of both the bosonic and fermionic 
condensates has always a zero at the origin. 

Next we study the stability of the solitons illustrated in figures 1 under a small 
perturbation inflicted by a sudden change in the boson-fermion scattering length cibf- 
After the solitons of figures 1 are formed we increase clbf by 10% at t = so that 
the nonlinearities Nbf and Nfb are suddenly increased from 40 to 44. This can be 
performed experimentally by varying a background magnetic field near a Feshbach 
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Figure 3. The stationary function |(^(,z,t)| for bosonic and fermionic dark solitons 
vs. z with Np — Nb — 100, v = 0.1 and asB = "bf = 5 nm in the presence 
of an optical-lattice potential with Ao = 5 and i>o = 5. The nonlinearities are 
N BB = 20, N BF = 40, ./Vps = 40, N FF = 1275. 

resonance in the boson-fermion system [39]. The solitons then execute small breathing 
oscillation around a mean position. The snapshot of the soliton profiles shown in figures 
3 (a) and (b) under this perturbation demonstrates their stability. 

Finally, we calculate the dark solitons in the boson-fermion mixture with the 
parameters of figure 1 (a) in the presence of an optical-lattice potential with v o = Ao = 5. 
We plot in figure 3 the function <p(z, t) for the dark soliton in this case. The presence of 
the optical-lattice potential creates modulations in both fermionic and bosonic functions. 
Because of the very strong repulsion, the fermions tend to occupy the whole available 
region in space and do not allow the formation of pronounced notches with hollow 
region inside at each optical-lattice site. Consequently, the modulation in the fermionic 
function is less pronounced than the modulation in the bosonic function. The bosonic 
function of figure 3 has pronounced notches compared to the smooth function in figure 
1 (a), whereas the fermionic function of figure 3 with small modulations is qualitatively 
more similar to the fermionic component of figure 1 (a). We also studied free expansion 
of the dark solitons in this case. The central notch in the fermionic dark soliton does 
not also expand in this case. No new interesting physics emerges and the results are not 
shown here. 
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4. Summary 

We use a coupled set of time-dependent mean-field-hydrodynamic equations for a 
trapped degenerate boson-fermion mixture to study the formation of a fermionic dark 
soliton (vibrational excitation) in a DFG as stationary states. We calculate the 
stationary functions with a notch at the center for fermionic dark soliton of the boson- 
fermion mixture. The existence of a central notch in the wave function for a dark soliton 
is typical of vibrational excitation. We perform numerical simulation of the dark solitons 
for a harmonic as well as a harmonic plus optical-lattice traps. The simulation is started 
with a time evolution of the mean-field equations with the eigenfunction of the lowest 
excited state of the linear oscillator problem setting all the nonlinearities to zero. The 
nonlinearities are introduced slowly during time evolution and the iteration continued 
until convergence is obtained. We demonstrate the stability of the dark soliton upon 
the application of a perturbation while the soliton executes small breathing oscillation. 

The present time-dependent formulation also permits us to study non-equilibrium 
free expansion of the coupled degenerate boson-fermion mixture. One way to observe 
the notch experimentally is to allow the dark soliton to undergo free expansion. We find 
that for a repulsive boson-fermion interaction, after a free expansion the notch in the 
fermion function for a relatively small fermion number of 100 does not increase in size 
significantly so as to be easily observable. 

In the present investigation we used a set of mean-field equations for the DFG- 
BEC mixture. A proper treatment of the DFG should be performed using a fully 
antisymmetrized many-body Slater determinant wave function [9]. Although we believe 
the present conclusion about fermionic dark soliton to be true in general, it would be of 
interest to establish their existence using such a fully antisymmetrized fermionic many- 
body wave function in the future. Also, a dark soliton features a particular spatial phase 
distribution, i.e. a step of n phase. In a BEC, this phase profile is supported due to the 
macroscopic phase of the condensate. Although a macroscopic phase of the fermionic 
component should emerge in the present mean-field model, to the best of our knowledge 
its existence in a fermionic many-body wave function has not been established rigorously. 
However, such a study is beyond the scope of this paper and would be a work of future 
interest. Nevertheless, it would be proper to call the fermionic excitation of the present 
paper by the term dark soliton due to its appropriate density distribution with a central 
notch. 
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